#coding:utf-8

from math import sqrt

xi = [7,13,15,25,31]
yi = [7,13,15,25,31]
zi = [7,13,20,25,31]
ri = [5,2,1,3,2]

n = 5
def soma_v(v1,v2):
	res = [0]*len(v1)
	for i in range(0,len(v1)):
		res[i] = v1[i]+v2[i]
	return res
	
def interno_v(v1,v2):
	res = 0
	for i in range(0,len(v1)):
		res += v1[i]*v2[i]
		
	return res	
	
def escalar_v(v1,k):
	res = [0]*len(v1)
	for i in range(0,len(v1)):
		res[i] = v1[i]*k
		
	return res	
	
def norma(v):
	res = 0
	for i in range(0,len(v)):
		res+=pow(v[i],2)
	res = sqrt(res)
	return res
	
def Ei(i,v,x,y,z):
	ind = 3*i
	return pow( (x-v[ind]),2 ) + pow( (y-v[ind+1]),2 ) + pow( (z-v[ind+2]),2 ) - pow( (ri[i]),2 )

def der_Ei_k(i,k,v,x,y,z):
	ind_i = 3*i	
	ind_k = 3*k
	
	dx = 2*(v[ind_k]-v[ind_i])
	dy = 2*(v[ind_k+1]-v[ind_i+1])
	dz = 2*(v[ind_k+2]-v[ind_i+2])
	
	return [ dx,dy,dz ]
	
def Fi(i,v):
	ind = 3*i
	res = 0
	for k in range(0,n):
		if k!=i:
			res+=Ei(k,v,v[ind],v[ind+1],v[ind+2])
	
	return res
	
def der_Fi_k(i,k,v):
	ind_i=3*i
	ind_k=3*k
	
	res = [0]*3
	
	if i==k:
		for j in range(0,n):
			if j!=i:
				soma_v(res,der_Ei_k(j,i,v,v[ind_i],v[ind_i+1],v[ind_i+2]))
	else:
		dx = -2*(v[ind_i]-v[ind_k])
		dy = -2*(v[ind_i+1]-v[ind_k+1])
		dz = -2*(v[ind_i+2]-v[ind_k+2])
		
		res[0] = dx
		res[1] = dy
		res[2] = dz
		
	return res	
	
def F(v):
	res = 0
	for i in range(0,n):
		res+=Fi(i,v)
		
	return res
	
def der_F_k(k,v):
	res = [0]*3
	for i in range(0,n):
		dfik = der_Fi_k(i,k,v)
		res[0]+=dfik[0]
		res[1]+=dfik[1]
		res[2]+=dfik[2]
		
	return res	
	
def grad_F(v):
	res = [0]*3*n
	for k in range(0,n):
		ind = 3*k
		dfk = der_F_k(k,v)
		res[ind] = dfk[0]
		res[ind+1] = dfk[1]
		res[ind+2] = dfk[2]
		
	return res
	
def v_from_vec():
	v = [0]*3*n
	for i in range(0,n):
		ind = 3*i
		v[ind] = xi[i]
		v[ind+1] = yi[i]
		v[ind+2] = zi[i]
		
	return v

def busca_linear(xk,dk):
	sigma = 0.8
	lbda = 1.0
	lim = 10
	
	while( F( soma_v(xk,escalar_v(dk,lbda)) ) >= F(xk) + sigma*lbda*interno_v(grad_F(xk),dk) ):
		if lim == 0:
			break			
		lbda/=2
		lim-=1
	
	return lbda
		
def main():
	x0 = v_from_vec()
	epson = 1e-2
	it = 0
	
	while(norma(x0)>epson):
		if it==100:
			break
			
		dk = escalar_v(grad_F(x0),-1)
		lbda = busca_linear(x0,dk)
		x0 = soma_v(x0,escalar_v(dk,lbda))
		
		it+=1
	
	for i in range(0,n):
		ind = 3*i
		print "P",i,"= ( ",x0[ind],x0[ind+1],x0[ind+2],ri[i]," )"
	
main()
